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A  NUMERICAL  SOLUTION  TECHNIQUE  FOR  MULTITEMPERATURE 
PLASMA  HYDRODYNAMICS 

INTRODUCTION 

There  are  many  problems  In  computational  plasma  dynamics  in  which  a 
multitemperature  fluid  treatment  of  the  plasma  is  sufficiently  complex  to 
adequately  describe  the  phenomena  of  primary  interest,  yet  is  sufficiently 
simple  to  permit  a  thorough  study  of  the  phenomena  without  a  prohibitive 
expenditure  of  resources.  In  a  multitemperature  approach  the  various 
particle  species,  and  perhaps  the  radiation  field,  are  modeled  as  distinct, 
interacting  fluids  which  exchange  energy  at  a  finite  rate  through  collisional 
or  collisionless  processes,  while  the  exchange  of  momentum  is  assumed  to  be 
sufficiently  rapid  that  the  velocities  of  the  component  fluids  are  identical. 
The  frequent  application  of  these  models  to  laboratory  and  astrophysical 
time-dependent  plasma  phenomena  make  valuable  a  computational  approach  to 
solving  the  equations  which  is  accurate,  efficient,  and  of  general  utility. 

We  describe  a  numerical  method  for  solving  the  multitemperature  plasma 

hydrodynamics  equations  in  a  timestep-split ,  finite-difference  formulation. 

Our  algorithm  possesses  several  important  characteristics.  We  solve  the 

convective  transport  equations  for  the  mass,  momentum,  and  total  energy 

densities  in  conservation  form,  thereby  guaranteeing  that  the  plasma  mass, 

momentum,  and  total  energy  be  locally  (and  therefore  globally)  conserved. 

At  the  same  time,  we  also  solve  the  convective  transport  equations  for  the 

internal  energy  densities  of  each  of  the  individual  species,  renormalizing 

the  resulting  solutions  so  that  they  are  consistent  with  the  total  energy 

equation.  In  this  way,  we  simultaneously  achieve  local  conservation  of  total 

energy  and  comparable  accuracy  in  the  description  of  each  of  the  plasma 
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species.  This  solution  of  an  overdetermined  set  of  hydrodynamic  equations  is 
a  unique  feature  of  the  algorithm.  Our  method  of  solution  yields  a  superior 
description  of  shock  wave  propagation  in  the  plasma,  also  as  a  consequence  of 
integrating  the  equations  in  conservation  form.  We  treat  the  temperature 
equilibration  between  species  in  a  fully  implicit  manner  in  the  diffusive 
transport  step,  and  as  a  result  equilibration  timescales  much  shorter  than 
the  hydrodynamic  timescales  in  the  problem  can  be  accomodated.  A  second 
novel  feature  of  the  algorithm  is  the  inclusion  of  the  temperature  changes 
induced  by  the  convection  as  source  terms  in  the  diffusive  transport 
calculation.  The  hydrodynamic  heating  then  acts  simultaneously  with  the 
other  processes  -  temperature  equilibration,  thermal  conduction,  etc.  -  which 
determine  the  species'  temperatures.  Finally,  for  problems  in  which  an 
implicit  calculation  of  the  thermal  conduction  terms  for  more  than  one  of  the 
plasma  species  is  necessary,  we  solve  the  coupled  temperature  equations  using 
an  efficient  iterative  scheme. 

The  simplest  example  of  a  multi  temperature  plasma  is  a  two-temperature 
model  of  a  plasma  composed  solely  of  electrons  and  ions.  This  model  contains 
all  of  the  essential  elements  of  the  general  approach,  and  we  therefore 
describe  in  detail  the  solution  procedure  for  this  simple  case.  We  conclude 
with  some  remarks  about  the  use  of  the  algorithm  and  its  generalization  to 
more  complex  plasma  models. 


THE  BASIC  EQUATIONS  OF  THE  TWO-TEMPERATURE  MODEL 


We  begin  by  considering  the  simple  case  of  a  two-component  plasma  of 
electrons  and  ions.  We  further  assume  that  the  plasma  is  charge-neutral  and 
that  the  exchange  of  momentum  between  the  electron  and  ion  fluids  is  perfect, 
so  that  the  plasma  is  also  current-free.  The  numerical  algorithm  that  we 
use  to  solve  this  problem  readily  generalizes  to  situations  involving  a 
larger  number  of  interacting  fluids.  The  four  basic  equations  of  the  model 
[1,2]  are  the  transport  equations  for  the  mass  density  p,  the  momentum 
density  pv,  the  electron  internal  energy  density  Ue>  and  the  ion  internal 
energy  density  U^, 


jf  +  V-(pv)  -  0,  (1) 

■—-(pv)  +  7  •  ( p vv )  +  7P  -  0,  (2) 

3U 

■t—  +  7»  (U  v)  +  P  ( V*v)  -  S  ,  (3) 

ot  e  e  e 

and 

3U 

+  V*  (U jV)  +  Pi(7*v)  =*  Si«  (4) 


The  other  quantities  in  these  equations  are  the  partial  pressures  Pg  and  P^ 
and  the  energy  source  terms  Sg  and  S^.  We  can  combine  Eqs.  (l)-(4)  to  obtain 

1  9 

a  fifth  equation,  for  the  total  energy  density  E  *  -^pv^  +  Ug  +  U^: 

f|  +  7* ( [E+P] v)  =  S,  (5) 
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where  P  ”  Pg  +  is  the  total  pressure  and  S  =  Sg  +  is  the  source  term 
for  the  total  energy.  Specification  of  equations  of  state  for  the  two 
species  completes  the  set  of  fundamental  relations  of  the  model. 

The  energy  source  terms  in  Eqs.  (3)-(5)  describe  the  local  deposition 

and  dissipation  of  energy  due  to  any  processes  relevant  to  the  physical 
problem  of  interest,  e.g.,  temperature  equilibration  between  species,  the 
conduction  of  heat  through  the  plasma,  bremmstrahlung  radiation  losses,  etc. 
The  thermal  coupling  of  the  two  fluids  leads  to  an  energy  exchange  rate  X,  so 

we  rewrite  the  electron  and  ion  source  terms  Sg  and  as 

S  -  Q  -  X  (6) 

e  e 

and 

Sj  -  Qi  +  X.  (7) 

We  adopt  for  the  exchange  rate  X  its  form  for  an  ideal  plasma  [1,2], 

X  -  C  u  (T  -  T  ) ,  (8) 

v  e  i 

in  which  Cv  is  the  specific  heat  capacity  at  constant  volume  for  an  ideal 
electron  gas,  w  is  the  reciprocal  of  the  temperature  equilibration  time,  and 
Tg  and  T^  are  the  electron  and  ion  temperatures. 
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THE  SOLUTION  TECHNIQUE  FOR  THE  TWO-TEMPERATURE  MODEL 

We  now  proceed  to  solve  the  initial-boundary  value  problem  posed  by  Eqs. 
(1)— (4),  together  with  Eqs.  ( 6 )— ( 8 )  for  the  energy  source  terms.  Our  method 
of  solution  is  based  on  an  Eulerian,  finite-difference  discretization  of  the 
equations,  in  which  we  timestep-split  the  convective  transport  from  the 
dissipative  processes.  The  splitting  permits  us  to  employ  an  explicit 
numerical  scheme  to  solve  the  convective  transport  problem  and  an  implicit 
scheme  to  solve  the  diffusive  transport  problem.  Significant  computational 
savings  can  result  from  this  approach  if  the  dissipative  processes, 
specifically  in  this  case  the  temperature  equilibration,  have  characteristic 
timescales  shorter  than  the  hydrodynamic  timescale  determined  by  the  Courant- 
Friedrichs-Lewy  stability  condition  [3],  If  the  convective  and  diffusive 
timescales  are  greatly  disparate,  however,  the  solutions  which  are  obtained 
should  be  tested  for  their  sensitivity  to  the  numerical  timestep  employed. 

In  the  convective  transport  step,  we  solve  the  overdetermined  set  of 
hydrodynamic  equations,  Eqs.  Cl)— (5),  setting  Sg  =  =  S  =  0.  We  use  a 

flux-corrected  transport  [4]  algorithm  to  integrate  the  equations.  The 
positivity  property  of  flux-corrected  transport  gives  it  an  advantage  over 
alternative  conservative  numerical  schemes,  e.g.  Lax-Wendroff  [5]  or 
conservative  leapfrog  [6],  which  may  in  principle  also  be  used.  Given  values 
for  all  of  the  physical  variables  at  time  t  (step  n),  we  begin  by  integrating 
these  equations  forward  to  time  t  +  At  (step  n+1),  symbolically 

r  n  ,  >.n  _n  ,,n  ,,ni 

Ip  »  (pv)  ,  E  ,  ug  ,  ul}  * 

5 
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E 


n+1 


i 

j  • 


r  n+1 

Ip 


/  +\n+l 
(pv) 


,  u 


u 


"n+1 
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We  may  now  calculate  the  total  internal  energy  density  U  either  as  the  sum  of 
the  partial  internal  energy  densities  Ug  and  ,  or  as  the  difference  between 
the  total  energy  density  and  the  kinetic  energy  density, 


,,'n+l  _  „'n+l  [(pv)n+1]2 

U  3  E  -  - — — . 

2pn+1 


(9) 


These  two  results  are  of  course  not  numerically  identical,  owing  to  the 
truncation  errors  in  the  differencing  of  the  convective  and  source  terms  in 
the  equations.  We  adopt  Eq.  (9)  as  the  correct  expression  for  the  total 
internal  energy  density,  thereby  ensuring  that  energy  is  conserved  locally  in 
the  convective  transport  calculation.  As  a  consequence,  we  must  specify  a 
prescription  for  correcting  the  partial  internal  energy  densities  so  that 
their  sum  is  equal  to  the  total  given  by  Eq.  (9). 


We  distribute  the  numerical  error  in  the  calculation  of  the  total 
internal  energy  density  locally  among  the  partial  internal  energy  densities 
in  the  following  way.  Assigning  fractional  errors  and  to  the  electrons 
and  ions,  we  calculate  at  each  mesh  point 


u’’n+1)], 


Requiring  that  ee  +  we  have 
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» 


J 

i 


as  desired.  A  prescription  for  distributing  the  error,  i.e.,  an  expression 

"n+l 

for  e  ,  is  to  be  selected.  Our  choice  is  to  take  e  proportional  to  U  ; 
s  s  s 

this  amounts  to  assuming  the  same  percentage  error  for  each  of  the  partial 
internal  energy  densities,  and  leads  to  the  renormalization  formulae 


'n+l 


U 


"n+l 


'n+l 


u"n+1  +  u"n+1 

e  i 


e,  i. 


(12) 


Alternative  prescriptions  for  distributing  the  error  are  certainly  available, 

e.g. ,  assuming  a  fixed  percentage  error  for  the  internal  energy  per  particle 

and  setting  e  proportional  to  the  local  number  density  n  .  We  prefer  the 
s  s 

prescription  above  because  of  the  simplicity  of  Eqs.  (12)  for  the  corrected 
partial  internal  energy  densities. 


By  solving  the  conservation  equations  for  the  mass,  momentum,  and  total 
energy  densities,  we  guarantee  local  conservation  of  mass,  momentum,  and 
energy  in  the  calculation,  and  we  ensure  that  the  Rankine-Hugoniot  conditions 
are  satisfied  across  any  shocks  that  develop  in  the  plasma.  At  the  same 
time,  we  calculate  the  internal  energy  densities  for  the  various  species 
independently  of  one  another,  so  that  the  state  of  each  of  the  plasma 
constituents  is  known  with  roughly  the  same  accuracy.  There  are  two  obvious 
alternatives  to  this  approach. 


The  first  alternative  is  to  solve  the  total  energy  equation  and  the 
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internal  energy  equation  for  one  of  the  species,  then  compute  the  internal 
energy  of  the  other  species  by  subtraction.  This  is  equivalent  to  setting 

I 

the  fractional  error  to  unity  for  the  second  species.  The  internal  energy 
density  obtained  in  this  manner,  and  therefore  the  pressure  and  temperature 
derived  from  it,  is  prone  to  large  errors  when  the  various  energies  in  the 

» 

problem  differ  by  orders  of  magnitude.  Our  algorithm  is  certainly  not  immune 
to  such  difficulties,  since  we  calculate  the  total  internal  energy  density  by 
subtraction  (Eq.  (9)).  Fortunately,  this  problem  is  self-correcting  as  the 
calculation  progresses,  owing  to  the  compressional  and  conduction  heating 
induced  by  the  abnormally  low  pressures  and  temperatures  in  the  regions  where 
the  internal  energy  becomes  negative.  Our  method  of  solution  at  least  has 
the  advantage  of  describing  each  of  the  species  with  comparable  accuracy 
throughout  the  plasma. 

» 

The  second  alternative  is  to  solve  only  the  equations  for  the  partial 
internal  energies.  This  approach  in  general  requires  specification  of  an 
artificial  viscosity  term  [7]  and  the  inclusion  of  viscous  stresses  in  the  I 

momentum  equation  and  viscous  heating  in  the  internal  energy  equation  for  at 
least  one  of  the  plasma  species.  By  judicious  choice  of  the  viscosity 

coefficients,  global  conservation  of  energy  can  be  achieved.  Our  approach  I 

has  the  advantage  of  not  requiring  the  introduction  of  an  artificial 

viscosity,  and  yet  guarantees  local  energy  conservation  in  a  straightforward 

manner.  Furthermore,  it  is  quite  generally  more  satisfactory  for  problems  in  » 

which  shock  waves  are  present.  Numerical  studies  [8,9]  show  that  the  shock 

speed  and  profile  are  consistently  more  accurately  determined  in  an  energy 

conservative  formulation  than  in  a  nonconservative  formulation.  » 
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We  now  turn  to  the  diffusive  transport  calculation.  In  the  diffusive 
transport  step  the  plasma  is  held  stationary  (v  =  0)  and  the  mass  density 
remains  at  its  value  at  the  end  of  the  convective  transport  step.  As  the 
intermediate  values  for  the  diffusive  transport  calculation  we  take  the  mass 
density  pn+*  at  the  end  of  the  current  convective  transport  step  and  the 
temperatures  T^11  and  T^n  from  the  preceding  diffusive  transport  step.  We 
formally  treat  the  temperature  changes  induced  by  the  convection  -  the 
hydrodynamic  heating  rates  -  as  additional  source  terms  in  the  energy 
aquations,  subject  to  redistribution  among  the  species  by  the  thermal 
coupling  between  the  fluids.  We  adopt  this  procedure  in  order  to  avoid 
potential  convergence  difficulties,  arising  in  regions  where  the  temperature 
equilibration  is  very  fast,  caused  by  allowing  the  convection  to  drive  the 
temperatures  apart  in  the  convective  step  and  then  having  the  temperature 
equilibration  attempt  to  force  them  back  together  in  the  subsequent  diffusive 
step. 


With  the  Inclusion  of  the  rates  of  hydrodynamic  heating  Hg  and  H^,  the 
equations  for  the  diffusive  transport  calculation,  obtained  from  Eqs.  (3)  and 
O),  are 


C  t —  =  Q  -  C  W  (T  -  T.)  +  H 

vedt  e  v  e  1  e 


and 
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"igure  2.  The  mass  densities  and  temperatures  inside  the  critical  surface  at 
the  peak  of  the  laser  pulse.  The  mass  densities  are  shown  in  (a), 

for  the  two-temperature  ( - )  and  one-temperature  ( - )  models. 

The  temperatures  are  shown  in  (b),  for  the  electrons  ( - )  and  the 

ions  (-•-)  in  the  two-temperature  model  and  for  the  one- 
temperature  model  ( - ). 
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Figure  l.  The  mass  densities  and  temperatures  in  the  laser-irradiated  target 
at  the  peak  of  the  laser  pulse.  The  mass  densities  are  shown  in 

(a),  for  the  two-temperature  ( - )  and  one-temperature  ( - ) 

models.  The  temperatures  are  shown  in  (b),  for  the  electrons 

( - )  and  the  ions  (-•-)  in  the  two-temperature  model  and  for  the 

one-temperature  model  ( - ). 
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individual  temperatures  (corresponding  to  Eqs.  (20)  and  (21))  are  solved 
iteratively  in  tandem  with  the  2(N— 1)  equations  for  the  difference 
temperatures  and  their  source  terras  (corresponding  to  Eqs.  (22)  and  (23)). 
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DISCUSSION 


The  numerical  algorithm  that  we  use  to  solve  the  equations  of  the  two- 
temperature  plasma  hydrodynamic  model  requires  in  the  convective  transport 
step  the  integration  of  the  equation  for  the  total  energy,  in  addition  to 
those  for  the  mass,  momentum,  and  partial  internal  energies.  The  cost  is  an 
increase  of  roughly  25%  in  the  amount  of  numerical  work  done  and  in  the 
computing  time  over  methods  which  neglect  either  the  total  energy  equation  or 
one  of  the  partial  internal  energy  equations.  The  benefits  are  a  guarantee 
of  local  conservation  of  total  energy,  a  more  accurate  description  of  the 
propagation  of  shock  waves  through  the  plasma,  and  comparable  accuracy  in  the 
computed  fluid  properties  of  each  of  the  plasma  species.  In  the  diffusive 
transport  step,  the  technique  of  iterating  between  equations  for  the  species' 
temperatures  and  equations  for  the  difference  temperatures  is  useful  for 
problems  in  which  the  conduction  of  heat  by  both  species  is  sufficiently 
rapid  to  call  for  an  implicit  treatment.  In  simpler  situations  a  direct 
solution  of  the  coupled  temperature  equations  is  possible. 

The  extension  of  the  algorithm  to  more  general  multiteraperature  plasma 
hydrodynamic  models  is  straightforward.  For  an  N-fluid  model,  the  convective 
transport  problem  consists  of  the  three  equations  for  the  mass,  momentum,  and 
total  energy  densities,  and  N  equations  for  the  partial  internal  energy 
densities.  After  calculating  the  total  internal  energy  density  from  the 
solution  of  the  first  three  equations,  the  separately  computed  partial 
internal  energy  densities  are  renormalized  by  an  appropriate  generalization 
of  Eqs.  (12).  In  the  diffusive  transport  step,  the  N  equations  for  the 
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plot  density  and  temperature  profiles  for  the  expanding,  accelerating  target 
at  the  peak  of  the  pulse  in  Figures  la,b.  Profiles  from  a  one-temperature 
model  calculation  are  shown  for  comparison  with  the  two-temperature  results. 
The  temperature  equilibration  time  varies  from  about  50  ns  in  the  blowoff 
plasma  to  less  than  1  ps  in  the  core.  The  ion  temperature  peaks  at  400  ev 
near  the  critical  surface,  where  thermal  coupling  between  the  electrons  and 
ions  is  still  effective,  and  then  declines  outward  from  that  point  as  the  ion 
fluid  cools  by  expansion.  The  electron  temperature  exhibits  a  monotonic 
increase  to  1100  ev  at  the  edge  of  the  foil,  due  to  inverse  bremmstrahlung 
absorption  in  the  underdense  plasma.  We  plot  the  density  and  temperature 
profiles  inside  the  critical  surface  in  Figures  2a, b.  The  electron  and  ion 
temperatures  converge  rapidly  toward  the  core  of  the  target,  coalescing  at 
about  100  ev  and  decreasing  together  to  less  than  1  ev  within  the  ablation 
surface. 
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A  LASER-MATTER  INTERACTION  SIMULATION 


We  have  implemented  the  two-temperature  algorithm  in  the  one-dimensional 
numerical  simulation  code  FASTID  [13].  This  code  was  developed  to  study  the 
plasma  hydrodynamics  of  the  interaction  of  intense  laser  radiation  with 
material  targets.  The  energy  in  the  incident  laser  beam  is  absorbed 
primarily  by  the  electrons  in  the  target,  which  then  equilibrate  locally  with 
the  ions  and  convect  and  conduct  energy  out  of  the  deposition  region. 
Conditions  in  the  irradiated  material  range  from  strong  thermal  coupling  in 
the  compressed,  shock-heated  core  of  the  target  to  negligible  coupling  in  the 
rapidly  expanding  blowoff  plasma.  In  the  simulations,  we  routinely  achieve 
mass,  momentum,  and  total  energy  conservation  to  within  machine  precision  and 
agreement  between  the  electron  and  ion  temperatures  in  the  strongly  coupled 
portion  of  the  target  to  better  than  one  part  in  one  thousand.  On  those 
occasions  when  the  total  internal  energy  in  a  computational  zone  is  driven 
negative  owing  to  numerical  errors,  the  algorithm  smoothly  corrects  the 
problem  over  several  timesteps  without  difficulty.  These  simulations  have 
tested  the  algorithm  over  an  extreme  range  of  physical  conditions.  Its 
robust  performance  suggests  that  it  is  a  reliable  method  of  solution  for 
general  hydrodynamic  problems  in  multitemperature  plasmas. 

We  simulated  the  response  of  a  6.5  urn  polystyrene  (CH)  planar  target  to 
a  gaussian  laser  pulse  of  wavelength  1.06  um,  peak  intensity  2  x  1013  W/cm2, 
and  full-width  half-maximum  duration  of  2.6  ns.  The  target  is  initially  at 
rest  at  x  *  0,  at  room  temperature  and  solid  density  (1.03  gra/cra3).  The  peak 
of  the  laser  pulse  is  reached  5  ns  after  the  start  of  the  calculation.  We 
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(26)  and  (27).  This  recalculation  is  necessary  in  order  to  maintain 
consistency  between  the  temperatures  used  in  the  diffusive  transport  step  and 
the  internal  energies  used  in  the  convective  transport  step.  The  integration 
of  the  plasma  hydrodynamic  equations  for  the  current  timestep  is  now 
complete. 


the  amount  of  computational  work  required,  this  conclusion  being  perhaps 
problem-  and  hardware-dependent. 


We  complete  the  diffusive  transport  calculation  by  updating  the  partial 
internal  and  total  energy  densities  and  the  temperatures.  Calculating  the 
accumulated  energy  changes  since  the  end  of  the  convective  transport  step. 


AU  ■  C  (Tn+lt“)  -  Tn)  -  H  At 
e  ve  e  e  e 


(24) 


and 


% "  cvi(Tr1(jo)  ■  ■  vc’ 


where  (°°)  denotes  the  value  at  the  end  of  the  temperature  iteration,  we 
calculate 


jn+1  =  U  n+l  +  AU  , 
e  e  e’ 


(25) 


if1  -  a;11*1  ♦ 


(26) 


(27) 


and 


En+1  -  £  n+l  +  AU  +  AU,. 

e  i 


(28) 


We  use  these  results  to  calculate  x  from  Eq.  (23), 

r  n+l(m)  n+l(m)i  n+l(m) 

J  X  y 

and  then  calculate  t  from  Eq.  (22), 

n+1 (m)  n+1 (m) 

X  T  • 

Finally,  substituting  the  expressions  for  and  into  Eqs.  (20)  and  (21), 

we  solve  them  separately  for  Tg  and  T^, 

Tn+l(m)  +  |Tn+l(m+l)  Tn+l(m+l)| 

The  iteration  cycle  may  now  be  carried  out  again,  starting  with  the  updated 
values  for  the  temperatures  Te  and  T^,  and  repeated  until  convergence  is 
achieved. 

As  an  alternative  to  the  iterative  method,  we  may  solve  the  coupled 
implicit  equations  for  the  temperatures  directly,  using  solution  techniques 
for  banded  or  sparse  matrices.  These  direct  methods  require  more  computer 
memory  and  time  to  solve  the  equations  than  does  the  iterative  method  to 
cycle  once,  but  the  direct  methods  are  of  course  one-pass  operations.  We 
have  found  in  practice  however  that  the  iteration  scheme  converges  very 
rapidly,  typically  to  a  tolerance  of  less  than  1%  in  two  or  three  cycles. 
Thus,  the  iterative  method  may  have  an  advantage  over  the  direct  methods  in 
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Dividing  Eqs.  (13)  and  (14)  for  the  species  temperatures  by  their 


respective  heat  capacities  and  subtracting  them,  we  obtain 


It  •  *  ~  aT« 


where  the  difference  temperature  t  is  defined  by 


t  -  Te  -  T1, 


the  source  term  x  is 


Qe+Hi  Qj  +  Hi 

C  C  .  ’ 

ve  vi 


and  the  relaxation  rate  0  is 


C  +  C  . 
ve  vi  _ 
ft  -  — - — -  C  (*>• 

C  C  .  v 
ve  vi 


We  difference  Eqs.  (13),  (14),  and  (16)  treating  the  difference  temperature  T 
fully  implicitly,  and  for  the  differencing  of  the  thermal  conduction  terras  in 


the  energy  sources  Qgand  Q^,  we  also  choose  a  fully  implicit  treatment  [11], 
Christiansen  and  Roberts  [10]  describe  a  time-centered  implicit  treatment  of 
the  equilibration  and  diffusion  [12]  processes.  Fixing  the  heat  capacities 


and  the  inverse  equilibration  time  at  their  values  at  the  end  of  the 
convective  transport  step,  we  obtain 


I 


where  Cvg  and  Cv^  are  the  specific  heat  capacities  at  constant  volume.  The 
hydrodynamic  heating  rates  are  given  in  finite-difference  form  as 


_  n+1  _  n 

,  , .  T  -  T 
,  n+1  s  _  s 


vs 


At 


e,  i. 


(15) 


in  which  we  calculate  the  heat  capacities  and  temperatures  at  step  n+1  from 

the  internal  energy  densities  U  and  U.  at  the  end  of  the  convective 

e  i 

transport  step.  We  calculate  the  intermediate  values  for  the  internal  and 


total  energy  densities  by  subtracting  the  energy  changes  HgAt  and  H^At 

*  *  * 

associated  with  the  compressional  heating  from  Ue>  U^,  and  E  (see  Eqs.  (24)- 
(28)  below). 


We  seek  solutions  of  Eqs.  (13)  and  (14)  in  which  the  temperature 
equilibration  is  treated  implicitly.  If  the  thermal  conduction  processes  may 
be  treated  explicitly,  then  we  trivially  solve  Eqs.  (13)  and  (14)  as  two 
equations  in  two  unknowns,  at  each  mesh  point.  If  an  implicit  treatment  of 
the  thermal  conduction  is  required  for  only  one  species,  we  can  eliminate  the 
temperature  of  the  other  species  algebraically  and  solve  the  resulting  system 
of  equations  by  standard  techniques  (e.g.,  tridiagonal  matrix  inversion  for  a 
three-point  diffusion  operator  in  one  space  dimension).  We  consider  the 
more  general  case  of  an  implicit  calculation  of  the  thermal  conduction  terms 
for  both  species,  to  which  we  apply  an  iterative  solution  technique.  Our 
approach  is  an  adaptation  of  the  iterative  method  described  by  Christiansen 
and  Roberts  [10]. 
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